Fast Monte Carlo algorithm for supercooled soft spheres 
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We present a nonlocal Monte Carlo algorithm with particle swaps that greatly accelerates thermal- 
ization of soft sphere binary mixtures in the glassy region. Our first results show that thermalization 
of systems of hundreds of particles is achievable, and find behavior compatible with a thermodynamic 
glass transition. 
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The dynamics and thermodynamics of supercooled liq- 
uids and glasses have received much attention for some 
time, but are still an open problem jjj. There has 
been theoretical progress regarding equilibrium dynam- 
ics H and thermodynamics (where older ideas |j|,Q have 
been recently developed PJq,|7|,p|,p|p0|) as well as non- 
equilibrium dynamics (mainly based on analogies with 
spin glass behavior p!l|). Many issues are open [^2[, such 
as the description of the dynamics of finite-range sys- 
tems, or the existence of a thermodynamic glass transi- 
tion, which is not accepted by all researchers [^3| . 

Given the great complexity of an analytical treatment, 
here as in other branches of science computer simulations 
play a major role |Q. However, simulations face difficul- 
ties, particularly in the study of equilibrium properties 
fl4j| . Well-known methods such as Molecular Dynamics 
(MD) and traditional Monte Carlo (MC), with dynamics 
that closely resembles that of real systems, suffer from a 
slow-down at temperatures approaching the glass tran- 
sition T g . Just as a glass well below T g remains out of 
equilibrium in the laboratory, model systems require pro- 
hibitively long CPU time to equilibrate, even for a hand- 
ful of particles. 

Thus in recent years several algorithms have been pro- 
posed endowed with a dynamics different from real sys- 
tems, one that would allow them to more easily jump 
across the free energy barriers that slow down MD or 
MC. We can mention replica MC p5| , multicanonical 
MC pq| , entropic sampling parallel tempering (PT) 
p8[ , expanded ensembles MC [j^), and the pivot clus- 
ter algorithm ||^,|2l[. Among these, parallel tempering 
is very convenient because of its straightforward imple- 
mentation and ease of parallelization, and because it can 
be applied to both MD and MC. It has been used suc- 
cessfully for spin glasses for several years, and recently 
applied to structural glasses p2j , p3| . However, even with 
these improvements, equilibration of deeply supercooled 
liquids requires a huge amount of CPU time for systems 
of N w 50 or more particles. 

In this paper we present a Monte Carlo algorithm with 
a simple but nonlocal dynamics which is suited for the 
study of the soft sphere binary mixture, a model struc- 
tural glass former widely studied in the past |24j2a,E0,E^| . 



This algorithm outperforms PT for soft spheres. It allows 
to equilibrate small systems with modest CPU time, and 
gives hope of equilibrating deeply supercooled systems of 
hundreds of particles. 

The soft spheres binary mixture is, with appropri- 
ately chosen parameters, a fluid in which crystallization 
is strongly inhibited [^4| . Half of the particles are of type 
A and half of type B, with radii a a and ob respectively 
and unit mass. The potential energy is 
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The radii are determined by fixing their ratio to be 
(jfl/o'/i = 1.2 and by setting the effective diameter to 
unity, which for equal number of A and B amounts to 
the condition 



(2a A ) 3 + 2(a a + a b f + (2a B ) 3 
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where Iq is the unit of length. Under these condi- 
tions, all excess thermodynamic properties depend only 
on T = pT _1/4 , with p = N/V the density and T the tem- 
perature. We have chosen Boltzmann's constant ks = 1, 
and we restrict ourselves to the case p = 1. The glass 
transition is known to happen at T c = 1.45 p5|| . Note 
that with this choices, the energy and temperature are 
dimensionless. 

We simulate this model with an algorithm that com- 
bines standard Monte Carlo moves with nonlocal moves 
(particle swaps) pq ], A step of the algorithm can be 
described as follows: 

• Choose a random particle i. 

• Draw 8 random numbers ri, . . . , r§ uniform in [0, 1]. 

• If r\ < 1 — ps then 

— Generate a standard MC trial move (i.e. shift 
particle by Ar = (2r 2 - 1, 2r 3 - 1, 2r 4 - l)Ar) 



else 



Choose a random particle j of a different type. 
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— Swap particle positions r$ <-> rj . 

— Shift both particles as in the standard move 
(Ari = (2r 2 - l,2r 3 - l,2r 4 - l)Ar, Ar, = 
(2r 5 -l,2r6-l,2r 7 -l)Ar). 

• Accept or reject new configuration according to 
the Metropolis criterion, i.e. accept it if rg < 
exp[-f3(E ncw - S oU )]. 

• Repeat for all particles. 

Ar and ps are the algorithm parameters. Ar is chosen 
at each temperature to achieve a shift acceptance ratio 
around 0.5, as is normally done in MC simulations, ps 
is the swap attempt probability, which we choose equal 
to 0.1. Implementation of this algorithm (which we shall 
call swap Monte Carlo (SW) algorithm) is straightfor- 
ward. Also, it should be obvious that SW satisfies de- 
tailed balance. 

Note that at high densities, even the swap of two neigh- 
boring particles is a cooperative (hence slow) process. 
The particles must go around each other, but they are 
hindered by the "cage" the rest of the particles build 
around them. The idea behind the SW algorithm is to ac- 
celerate these processes by making them non-cooperative. 

We have run SW for systems of N = 34 and N = 800 
particles in a cubic box with periodic boundary condi- 
tions at several values of T up to T = 2, deep in the 
glassy region. For the larger system, a long-range cutoff 
at r 2 = 3 was imposed. 

We start by comparing the approach to equilibrium of 
a system of N = 34 particles simulated with standard 
MC, parallel tempering MC and the SW algorithm. The 
PT data are from ref. [^|| and were obtained in a PT 
run with 13 replicas with T = 1, 1.05, . . . , 1.2, 1.3, . . . , 2. 
The evolution of the (potential) energy per particle for 
each algorithm starting from a random configuration is 
shown in figs. |l| and |^ for temperatures above and below 
the transition. At T = 1.4 (above the transition), all 
algorithms equilibrate rapidly. At T = 1.8, however, it 
is clear that relaxation times have soared, and the three 
algorithms behave differently. SW is seen to relax much 
faster than MC and PT. 
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FIG. 1. Comparison of the evolution the energy start- 
ing from a random configuration for the three algorithms at 
T = 1.4. 
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FIG. 2. As fig. but for V = 1.8. 

It may seem surprising that SW works at all, i.e. that 
the swaps are ever accepted. But the radii of the spheres 
are not very different from each other, and both species 
are present in equal amounts, so the replacement of one 
species by the other is not always very unfavorable. In- 
deed, the swap acceptance ratio is very small (at T = 1.8 
it is about 10 -4 ). However, this means that a small frac- 
tion of particles can escape the cage effect, with dramatic 
consequences for the dynamics. 

To check whether SW is capable of equilibrating the 
system we run simulations with 2 different initial config- 
urations at r = 1.4, 1.5, . . . , 2. We find that after 1-2 mil- 
lion steps (about 80 minutes of CPU time on an Alpha 
workstation), both configurations stabilize at the same 
energy. The runs are then continued up to 10 7 steps, 
using the last 8 • 10 6 steps to compute the energy and 
specific heat. We check thermalization in two ways: 

• We divide the last part of the run in blocks of 1 
million steps and verify that there is no systematic 
shift in the average energy of these blocks. 

• We compute the specific heat in two different ways, 
through C v = d(E)/dT and through C v = ((E 2 )- 
(E) 2 )/T 2 , and find agreement between the two. 

The energy and specific heat are shown in figs. || and |] 
as a function of temperature. The specific heat obtained 
from fluctuations agrees with the value from the deriva- 
tive of the energy, showing that the system has thermal- 
ized. A broad but clear peak around T w 0.12 is found 
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FIG. 3. Average energy vs. temperature for N = 34. Lines 
are only to join neighboring points. 
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FIG. 5. Energy vs. Monte Carlo time for N = 800, av- 
eraged over 4 initial configurations. From top to bottom, 
r = 1.4, 1.5, 1.6, 1.7, 1.8, 2. Lines are power law fits. 



> 
U 



1.9 



1.6 



1.5 





from fluctuations 
from derivative 


X 
X 


X 
X 


X 

X 

X 




X 

- / 


X 


X 


\/ 1 







0.05 



0.1 



0.15 0.2 
Temperature 



0.25 



0.3 



FIG. 4. Specific heat vs. temperature for N = 34. The 
values obtained from the derivative of the energy and from 
energy fluctuations fall on the same curve. 
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FIG. 6. Exponent a of the energy relaxation vs. tempera- 
ture. 



Next we simulate a system of TV = 800 particles at 
r = 1.2, 1.3, . . . , 1.8, 2. Fig. H show the energy relaxation 
for runs of 10 6 steps. The lowest temperatures clearly 
have not reached equilibrium in this case. The energy 
relaxation (after an initial fast decay which corresponds 
to the equilibration time for the liquid phase, namely 10 4 
steps) can be described by a power law E — E + t~ a , 
with a strongly temperature dependent (fig. ^J). A sharp 
change in slope of the a vs. T curve is observed near 
T w 0.15. The asymptotic energy (fig. also shows 
singular behavior at this temperature, as evidenced by 
the specific heat (fig. ^) estimated from the derivarive of 
the asymptotic energy. 
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FIG. 7. Asymptotic energy vs. temperature. Lines are only 
to join neighboring points. 
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FIG. 8. Specific heat vs. temperature for TV = 800 obtained 
from the asymptotic value of the energy (points) with theo- 
retical prediction from ref. J§| (continuous curve). 

Equilibrium simulations at low temperatures are rele- 
vant to the question of the existence of a thermodynamic 
glass transition driven by an entropy crisis p|,|3C|] . In this 
respect it is to note that the peak in the specific heat of 
the N = 34 system (fig. |) at T w 0.12 corresponds to 
T « 1.7, i.e. close to the predicted J§| critical temperature 
Tk ~ 1-65 at which a thermodynamic transition would 
take place and where the calculated specific heat shows 
a sharp drop to a value of 3/2 in the low temperature 
phase. 

Remarkably, a similar but sharper behaviour is found 
for the N = 800 system, though admittedly these values, 
resulting from extrapolation, should be taken with care. 
In figure || we also show the theoretical prediction for 
the specific heat (from Q). There is a difference in the 
liquid phase, which is to be expected since the hyper- 
netted chain approximation (HNC) used is known not to 
perform very well at low temperatures. But the position 
of the jump, at T « 1.7, is in good agreement with the 
theoretical Tk ~ 1-65 and with the position of the peak 
for the N = 34 system. The sharp change in the slope of 
the dynamic exponent a vs. T curve takes place at this 
same temperature. 

In summary, we have presented a very efficient MC 
algorithm for equilibration of a soft spheres binary mix- 
ture. We have managed to equilibrate an N — 34 system 
up to r = 2 in roughly 2 • 10 6 steps (about 80 minutes 
of CPU time on an Alpha workstation). Thus the SW 
algorithm puts thermalization of undercooled systems of 
hundreds of particles within reach of modern CPUs. Our 
first results with SW find behavior compatible with a 
thermodynamic glass transition. 

It would be very interesting to apply SW to other bi- 
nary systems. However, it is likely that it will need some 
modification to achieve a reasonable swap acceptance ra- 
tio. 
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